#######################################################################
#                                                                     #
#      Enhanced Creativity in Autism is Due to Co-occurring ADHD      #
#                                                                     #
#                   Emily C. Taylor, Gosia Goclowska,                 #
#              Mitchell J. Callan, & Lucy A. Livingston               #
#                                                                     #
#        Journal of Psychopathology and Clinical Science, 2024        #
#                                                                     #
#                         Script for use in R                         #
#######################################################################

library(psych) # Summaries by group
library(sjPlot) # Regression and t-test tables including bootstrapping
library(rstatix) # Cohens D with CIs
library(dplyr) # Filtering by AQ score
library(robustbase) # Robust regressions

data <- read.csv("Supplementary Data.csv", fileEncoding = "UTF-8-BOM")

####Main Text
###Table 1 Matching Autistic and Non-Autistic Groups
table(data$Sex, data$ASD)
describeBy(data[c("Age","ICAR", "AQ", "ASRS")], group=data$ASD)

t1 <- lm(Age ~ ASD, data = data)
t2 <- lm(ICAR ~ ASD, data = data)
t3 <- lm(AQ ~ ASD, data = data)
t4 <- lm(ASRS ~ ASD, data = data)

#Standard t-tests
tab_model(t1, t2, t3, t4, show.stat = TRUE, show.est = FALSE, string.stat = "t", show.intercept = FALSE, show.obs = FALSE, show.r2 = FALSE)

#Cohen's Ds
set.seed(1234)
cohens_d(Age ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(ICAR ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(AQ ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(ASRS ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)

###Table 2 Group Mean and Mean Differences in Creativity
describeBy(data[c("DT_Fluency","DT_Flexibility", "Sqrt_DT_Originality", "Log_CAQ", "BICB", "CPS", "CSE")], group=data$ASD)

t5 <- lm(DT_Fluency ~ ASD, data = data)
t6 <- lm(DT_Flexibility ~ ASD, data = data)
t7 <- lm(Sqrt_DT_Originality ~ ASD, data = data)
t8 <- lm(Log_CAQ ~ ASD, data = data)
t9 <- lm(BICB ~ ASD, data = data)
t10 <- lm(CPS ~ ASD, data = data)
t11 <- lm(CSE ~ ASD, data = data)

#Standard t-tests
tab_model(t5,t6,t7,t8,t9,t10,t11, show.stat = TRUE, show.est = FALSE, string.stat = "t", show.intercept = FALSE, show.obs = FALSE, show.r2 = FALSE)

#Bootstrapped t-tests
tab_model(t5,t6,t7,t8,t9,t10,t11, bootstrap = TRUE, seed = 1234, iterations = 5000, show.ci = 0.95, show.intercept = FALSE, show.obs = FALSE, show.r2 = FALSE, show.p = FALSE)

#Cohens Ds
set.seed(1234)
cohens_d(DT_Fluency ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(DT_Flexibility ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(Sqrt_DT_Originality ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(Log_CAQ ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(BICB ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(CPS ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(CSE ~ ASD, data = data, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)

### Table 3 Multiple Regression Analyses Exploring Autism and ADHD as Predictors of Creativity
##Analysis 1: Clinical Diagnoses
m1 = lm(DT_Fluency ~ ASD + ADHD + ICAR + Sex + Age, data = data)
m2 = lm(DT_Flexibility ~ ASD + ADHD + ICAR + Sex + Age, data = data)
m3 = lm(Sqrt_DT_Originality ~ ASD + ADHD + ICAR + Sex + Age, data = data)
m4 = lm(Log_CAQ ~ ASD + ADHD + ICAR + Sex + Age, data = data)
m5 = lm(BICB ~ ASD + ADHD + ICAR + Sex + Age, data = data)
m6 = lm(CPS ~ ASD + ADHD + ICAR + Sex + Age, data = data)
m7 = lm(CSE ~ ASD + ADHD + ICAR + Sex + Age, data = data)

#OLS Regressions
tab_model(m1,m2,m3,m4,m5,m6,m7, show.intercept = FALSE, show.se = TRUE, show.std = TRUE, show.obs = FALSE, show.r2 = FALSE, show.ci = FALSE)

#Bootstrapped Regressions
tab_model(m1,m2,m3,m4,m5,m6,m7, bootstrap = TRUE, seed = 1234, iterations = 5000, show.ci = 0.95, show.intercept = FALSE, show.se = TRUE, show.obs = FALSE, show.r2 = FALSE, show.p = FALSE)

##Analysis 2: Traits
m8 = lm(DT_Fluency ~ AQ + ASRS + ICAR + Sex + Age, data = data)
m9 = lm(DT_Flexibility ~ AQ + ASRS + ICAR + Sex + Age, data = data)
m10 = lm(Sqrt_DT_Originality ~ AQ + ASRS + ICAR + Sex + Age, data = data)
m11 = lm(Log_CAQ ~ AQ + ASRS + ICAR + Sex + Age, data = data)
m12 = lm(BICB ~ AQ + ASRS + ICAR + Sex + Age, data = data)
m13 = lm(CPS ~ AQ + ASRS + ICAR + Sex + Age, data = data)
m14 = lm(CSE ~ AQ + ASRS + ICAR + Sex + Age, data = data)

#OLS Regressions
tab_model(m8,m9,m10,m11,m12,m13,m14, show.intercept = FALSE, show.se = TRUE, show.std = TRUE, show.obs = FALSE, show.r2 = FALSE, show.ci = FALSE)

#Bootstrapped Regressions
tab_model(m8,m9,m10,m11,m12,m13,m14, bootstrap = TRUE, seed = 1234, iterations = 5000, show.ci = 0.95, show.intercept = FALSE, show.se = TRUE, show.obs = FALSE, show.r2 = FALSE, show.p = FALSE)

####Supplementary Material
###Supplementary Table 1 Across Sample Mean Scores and Distributions for Each Measure of Creativity
describe(data[c("DT_Fluency","DT_Flexibility", "DT_Originality", "CAQ", "BICB", "CPS", "CSE")])

###Supplementary Table 2 Autistic and Non-Autistic Group Demographics after applying Inclusion Criteria based on Autism Spectrum Quotient Thresholds
data_AQinclusion <- filter(data, (ASD == 1 & AQ >= 32) | (ASD == 0 & AQ < 26))

table(data_AQinclusion$Sex, data_AQinclusion$ASD)
describeBy(data_AQinclusion[c("Age","ICAR", "AQ", "ASRS")], group=data_AQinclusion$ASD)

t12 <- lm(Age ~ ASD, data = data_AQinclusion)
t13 <- lm(ICAR ~ ASD, data = data_AQinclusion)
t14 <- lm(AQ ~ ASD, data = data_AQinclusion)
t15 <- lm(ASRS ~ ASD, data = data_AQinclusion)

#Standard t-tests
tab_model(t12, t13, t14, t15, show.stat = TRUE, show.est = FALSE, string.stat = "t", show.intercept = FALSE, show.obs = FALSE, show.r2 = FALSE)

#Cohen's Ds
set.seed(1234)
cohens_d(Age ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(ICAR ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(AQ ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(ASRS ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)

###Supplementary Table 3 Group Mean and Mean Differences in Creativity after applying Inclusion Criteria based on Autism Spectrum Quotient Thresholds.

describeBy(data_AQinclusion[c("DT_Fluency","DT_Flexibility", "Sqrt_DT_Originality", "Log_CAQ", "BICB", "CPS", "CSE")], group=data_AQinclusion$ASD)

t16 <- lm(DT_Fluency ~ ASD, data = data_AQinclusion)
t17 <- lm(DT_Flexibility ~ ASD, data = data_AQinclusion)
t18 <- lm(Sqrt_DT_Originality ~ ASD, data = data_AQinclusion)
t19 <- lm(Log_CAQ ~ ASD, data = data_AQinclusion)
t20 <- lm(BICB ~ ASD, data = data_AQinclusion)
t21 <- lm(CPS ~ ASD, data = data_AQinclusion)
t22 <- lm(CSE ~ ASD, data = data_AQinclusion)

#Standard t-tests
tab_model(t16,t17,t18,t19,t20,t21,t22, show.stat = TRUE, show.est = FALSE, string.stat = "t", show.intercept = FALSE, show.obs = FALSE, show.r2 = FALSE)

#Bootstrapped t-tests
tab_model(t16,t17,t18,t19,t20,t21,t22, bootstrap = TRUE, seed = 1234, iterations = 5000, show.ci = 0.95, show.intercept = FALSE, show.obs = FALSE, show.r2 = FALSE, show.p = FALSE)

#Cohens Ds
set.seed(1234)
cohens_d(DT_Fluency ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(DT_Flexibility ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(Sqrt_DT_Originality ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(Log_CAQ ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(BICB ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(CPS ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)
set.seed(1234)
cohens_d(CSE ~ ASD, data = data_AQinclusion, ref.group = 2, var.equal = TRUE, ci = TRUE, nboot = 5000)

###Supplementary Table 4 Multiple Regression Analyses Exploring Autism and ADHD as Predictors of Creativity after applying Inclusion Criteria based on Autism Spectrum Quotient Thresholds

##Analysis 1: Clinical Diagnoses
m15 = lm(DT_Fluency ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)
m16 = lm(DT_Flexibility ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)
m17 = lm(Sqrt_DT_Originality ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)
m18 = lm(Log_CAQ ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)
m19 = lm(BICB ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)
m20 = lm(CPS ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)
m21 = lm(CSE ~ ASD + ADHD + ICAR + Sex + Age, data = data_AQinclusion)

#OLS Regressions
tab_model(m15,m16,m17,m18,m19,m20,m21, show.intercept = FALSE, show.se = TRUE, show.std = TRUE, show.obs = FALSE, show.r2 = FALSE, show.ci = FALSE)

#Bootstrapped Regressions
tab_model(m15,m16,m17,m18,m19,m20,m21, bootstrap = TRUE, seed = 1234, iterations = 5000, show.ci = 0.95, show.intercept = FALSE, show.se = TRUE, show.obs = FALSE, show.r2 = FALSE, show.p = FALSE)

##Analysis 2: Traits
m22 = lm(DT_Fluency ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)
m23 = lm(DT_Flexibility ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)
m24 = lm(Sqrt_DT_Originality ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)
m25 = lm(Log_CAQ ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)
m26 = lm(BICB ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)
m27 = lm(CPS ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)
m28 = lm(CSE ~ AQ + ASRS + ICAR + Sex + Age, data = data_AQinclusion)

#OLS Regressions
tab_model(m22,m23,m24,m25,m26,m27,m28, show.intercept = FALSE, show.se = TRUE, show.std = TRUE, show.obs = FALSE, show.r2 = FALSE, show.ci = FALSE)

#Bootstrapped Regressions
tab_model(m22,m23,m24,m25,m26,m27,m28, bootstrap = TRUE, seed = 1234, iterations = 5000, show.ci = 0.95, show.intercept = FALSE, show.se = TRUE, show.obs = FALSE, show.r2 = FALSE, show.p = FALSE)

###Supplementary Table 6 Robust Regression Analyses Exporing Autims and ADHD as Predictors of Creativity
##Analysis 1: Clinical Diagnoses
robm1 <- lmrob(DT_Fluency ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")
robm2 <- lmrob(DT_Flexibility ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")
robm3 <- lmrob(Sqrt_DT_Originality ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")
robm4 <- lmrob(Log_CAQ ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")
robm5 <- lmrob(BICB ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")
robm6 <- lmrob(CPS ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")
robm7 <- lmrob(CSE ~ ASD + ADHD + ICAR + Sex + Age, data = data, setting="KS2014")

#Regressions
tab_model(robm1,robm2,robm3,robm4,robm5,robm6,robm7, show.intercept = FALSE, show.se = TRUE, show.std = TRUE, show.obs = FALSE, show.r2 = FALSE, show.ci = 0.95)

##Analysis 2: Traits
robm8 <- lmrob(DT_Fluency ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")
robm9 <- lmrob(DT_Flexibility ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")
robm10 <- lmrob(Sqrt_DT_Originality ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")
robm11 <- lmrob(Log_CAQ ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")
robm12 <- lmrob(BICB ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")
robm13 <- lmrob(CPS ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")
robm14 <- lmrob(CSE ~ AQ + ASRS + ICAR + Sex + Age, data = data, setting="KS2014")

#Regressions
tab_model(robm8,robm9,robm10,robm11,robm12,robm13,robm14, show.intercept = FALSE, show.se = TRUE, show.std = TRUE, show.obs = FALSE, show.r2 = FALSE, show.ci = 0.95)